Nonlinear lattice model of viscoelastic Mode III fracture 
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We study the effect of generai nonlinear force laws in viscoelastic lattice models of fracture, 
focusing on the existence and stability of steady-state Mode III cracks. We show that the hysteretic 
behavior at small driving is very sensitive to the smoothness of the force law. At large driving, 
we find a Hopf bifurcation to a straight crack whose velocity is periodic in time. The frequency of 
the unstable bifurcating mode depends on the smoothness of the potential, but is very close to an 
exact period-doubling instability. Slightly above the onset of the instability, the system settles into 
a exactly period-doubled state, presumably connected to the aforementioned bifurcation structure. 
We explicitly solve for this new state and map out its velocity-driving relation. 



I. INTRODUCTION 



The problem of dynamic fracture has received increasing attention in the physics community in the last decade [El . 
The experiments of Fineberg, et. al. showing interesting dynamical behavior for large velocity cracks have been 
at the center of this growing interest. From a theoretical point of view, the singularities at the crack tip associated 
with continuum treatments of the crack problem make the problem challenging. The presence of these singularities 
necessitates a treatment of the crack at the microscopic level, where the continuum treatment and its associated 
singularities are not applicable. One line of attack on this problem, initiated by Slepyan ||, has been through the 
Q ■ study of lattice models of cracks. These lattice models are simpler than full atomistic simulations in that the 
O \ connectivity of the atoms is specified from the beginning, and so dislocations are excluded. 

Much progress has been made in understanding steady-state propagation of cracks in lattice systems [[§ |], @, f§ • 
I ' A key simplifying assumption underlying much of this progress has been the assumption of piecewise-linear forces 
, between the particles, so that the particles interact with Hookean springs, which break at some critical extension, 
0^ ' reducing the force to zero. With these piecewise-linear interactions, the model admits an analytic solution via the 
Wiener-Hopf technique. This solution has been carried out both for Mode I and Mode III cracks, for both finite width 
and infinitely wide systems, with and without dissipation. 

A general feature of these solutions is that they are problematic at both very small and very large (roughly above 
0.7 of the wave speed) velocities. For small dissipation, the solutions are inconsistent at small velocities, in that 
bonds on the crack surface which are assumed to crack at time t = 0, say, in fact are seen to crack earlier. For large 
dissipation, the small velocity solutions are consistent, but exhibit a backward (velocity increasing with decreasing 
driving) velocity/driving curve, indicating the solutions are unstable. Thus, in any case, stable solutions are not found 
for small velocities. 

At large velocities, the analytic solutions are again inconsistent, this time due to the breaking of additional bonds 
not on the crack surface. The analytic methods thus are unable to tell us anything about the dynamics beyond this 
point. 

An additional limitation of the piecewise-linear force law is that it complicates the task of constructing a linear 
stability treatment of the steady-state crack. This is due of course to the discontinuous nature of the force law. For 
these reasons, we choose in this paper to examine steady-state crack propagation in lattice models with arbitrary 
force-laws. In particular, we study a family of force-laws parameterized by a, such that for small a the force-law is 
smooth and the force-law goes over to the piecewise-linear one in the limit of infinite a. We have previously studied (^] 
the behavior of arrested Mode III cracks with this family of force-laws, and found that the range of drivings for which 
arrested cracks exists narrows sharply as a is reduced from infinity. In this paper we examine moving Mode III cracks 
for varying a. We find that the behaviour for small velocity is very sensitive to a. For large velocities, we find that 
the effect of finite a is to convert the inconsistency to a regular linear instability, here of Hopf type. We show that 
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at some distance into the unstable regime, the system adopts a new form of steady-state behavior which breaks the 
symmetry across the crack surface and has a period-2 structure. We conclude with some observations about the model 
deep in the unstable regime. 



II. MODEL AND GENERAL METHODOLOGY 



In this paper, we study a triangular lattice model of Mode III cracking. We take the force exerted at site X\ by the 
displacement at nearest neighbor site X2, with (scalar) displacements Ui, 112, respectively to be 

1 + tanh(a(l - cri. 2 (u 2 - wi))) 

/l,2 = (U2 ~ Ul) — — r—^ (1) 

1 + tanh(a) 

where a is a parameter which controls the smoothness of the potential and a is positive if X2 lies to the left of x\ or is 
further from the crack plane y = 0. This form has the feature that the force goes to zero at large positive extensions 
of the nearest-neighbor springs. In the limit of large a, this force-law approaches a step-function, (with critical spring 
extension 1), a form introduced by SlepyanQ, and the subject of much investigation^ ||, fj], |], [h], The lattice has 
2N + 2 rows in the y-direction, separated by a distance V3/2, so that the rows are labeled from y = -(2N + l)\/3/4 
to y — (2N + l)\/3/4. The displacements u of the bottom and top rows at y = ±(2JV + l)V3/4 are constrained to be 
±A. We introduce a Kelvin-type viscosity parameterized by r\ via an additional dissipative force 

01,2 = r]kx > 2^-{u2 - ui) (2) 
at 



where is an effective spring constant 



fel,2 = /l,2/(w2 - Ul) (3) 



This form was chosen so as to ensure a purely dissipative force which in the limit of large a goes over to the form 
r)9(l — (u2 — U\)){u2 — iii) studied in connection with dissipation in the piecewise- linear model §, @, 
These forces define the model. The equation of motion is simply 

Mg= ^ (/x,a» + 9x,x>) (4) 
x' £nn 

We shall primarily be interested here in steady-state cracks, where the displacements have the Slepyan form u XtV (t) = 
u y (t — x/v). Furthermore, we will focus on symmetric cracks, such that U—y(t) = —u v (t + l/(2w)), which is the 
symmetry appropriate to steady-state cracks in the piecewise-linear model B. Then, solving for the steady-state 
crack means solving for N functions of time, characterizing the time development of a typical mass point in each row 
of the lattice. The equations are non-local in time, due to the coupling between different lattice points. As there is 
no hope of constructing analytic solutions at general a, we solve this problem numerically. 

To proceed, we discretize the time variable in units of some small dt. The key insight involved in constructing a 
numerical procedure is noting that the steady-state equations have a banded structure, just as was the case for the 
piecewise-linear model Q. This is due to the fact that mass points are coupled to nearest neighbors, which gives rise 
via the Slepyan ansatz to a coupling to displacements at a finite time separation. Of course, the equations here are 
everywhere nonlinear and we need to employ Newton's method to find a solution. This is nonetheless computationally 
feasible, since the update step in Newton's method is a linear problem with the aforementioned banded structure. 

Imagine searching for a solution assumed to have some velocity v. There is one equation of motion for each 
displacement field u y at each discrete instant in time. In addition, fixing the translation invariance by some condition 
such as w^/g; 4 (0) = 1/2 gives one additional equation, which we can use to solve for the driving displacement A. The 
only problem is that A is involved in many different equations and so destroys the banded structure of the problem. 
Two different ways to circumvent this difficulty present themselves. The simpler of the two is to guess, for a given 
velocity, a value of A, and solve a modified system of equations where the equation of motion at t = 0, y = \/3/4 has 
been dropped, and the time translation has been fixed. The resulting banded problem can be solved economically, 
and the violation of the dropped equation of motion calculated. This defines a mismatch function, which has a zero at 
the true value of A, which is located by a standard zero-finding routine. A more efficient approach is to realize that 
solving the entire linear system can be accomplished at essentially no more expense than the fully banded modified 
system of the first approach. The algorithm for achieving this is presented in Appendix A. The advantage of this 
approach is that the Newton solver produces a value of A directly, without the need of invoking a outer root-finding 
routine. 
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III. THE SMALL VELOCITY REGIME 



We have already seen in that smoothing the form of the force-law effects a dramatic reduction in the window of 
drivings for which arrested cracks exist. As the moving crack solution arises as a backward bifurcation of the arrested 
crack, we can expect that small velocity cracks are also extremely sensitive to the smoothness of the force-law. An 
indication of this can be found in our previous study || of a continuous- a;, discrete-?/ model, which did not have any 
window of arrested cracks. There the velocity rose linearly from zero as the driving A increased from the Griffith value 
Ag, with a slope inversely proportional to 77. The study of small velocity solutions with sharp but not discontinuous 
force laws should be directly related to experimental data on the onset of propagating cracks in materials such as 
single crystal silicon |l3| . 

As mentioned above, we solve the steady-state problem via Newton's method. We present in Fig. 0(a-c) a graph of 
v versus A/Ag for various values of rj and a, together with the results from the piecewise-linear limit. We chose to 
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FIG. 1: Dependence of v on A/Aq for a = 5, 15, 40, and 100 together with the piecewise-linear model for 77 = 0.25, 1 and 3. 
N = 3, dt — 0.1. 



present data for N — 3 so as to be able to investigate smaller velocities at reasonable computational cost; the small 
velocity regime is relatively insensitive to N ||. A few general trends are evident from these plots. First, the effect 
of a is much more pronounced at small velocity. Second, smoothing the potential by decreasing a postpones the 
onset of the backward branch of the curve to lower velocity. This is consistent with the narrowing of the window of 
arrested cracks with decreasing a, since the curve is forced to turn back to meet the end of the arrested crack branch. 
Decreasing a also decreases the amplitude of the oscillations present at small 77. Examining the effect of varying 77, we 
see that the larger 77's are less sensitive to a. This is to be expected, since even in the piecewise-linear model r\ reduces 
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FIG. 2: Dependence of ryu on A/A G for a = 5, for r? = 0.25, 1, and 3. Again, N = 3, dt = 0.1. 



the extent of the backward branch which is one of the primary consequences of having smaller a. This reduction is 
related to the increase in the size of the process zone with increasing r\ , which smoothes out the lattice structure 
in a similar manner to that accomplished by smoothing the potential. Thus, in general, decreasing a and increasing 
77 both act to reduce the lattice-trapping effects which give rise to the backward branch. It should be noted that this 
analysis is consistent with the onset of running cracks directly above Aq in a molecular dynamics simulation using 
Lenard- Jones potentials [fu||. 

Another way to see that reasonably smooth force-laws give rise to essentially ai-continuum behavior is to look at 
-qv as a function of A for different 77's. In Fig. [|, we present the data for the case a = 5. We see that the data 
overlap extremely well, except for the highest velocity data for each value of 77. This is an indication that, as in 
the a;-continuum problem ||, the only velocity scales are rj and the wave speed, c = -\/3/2. For the piecewise-linear 
model, on the other hand, where the x-lattice scale is important, the rjv scaling only sets in for large 77, where 
the process zone is large 0. 

The other point of interest is the nature of the solutions for small 77. In the piecewise-linear model, the solutions 
on the backward branch are inconsistent at small 77, since the underdamping of the forward running waves leads to 
pre-cracking. This is of course not a problem in our fully nonlinear model, since cracking here is a reversible process. 
Nevertheless, it is amusing to note that the solutions on the backward branch for small 77 do not behave in this 
manner. They do not crack and reheal. Rather, they exploit the smooth transition in the force law to crack in a slow, 
monotonic fashion. We show in Fig. |^a the extension of the bond along the crack surface as a function of time, for 
a — 40, 100, and 200, along with the analogous result for the piecewise-linear model. We see that the limit of large 
a does not correspond to the behavior of the piecewise-linear model. Not surprisingly, then, it turns out the v vs. A 
curve for large a does not converge to the piecewise-linear result. Thus, for example, for the case considered in Fig. 
0a, namely N=3, v = .1, T] = 0.25, the limiting value of A/Aq for large a is 1.2140, as compared to 1.2778 for the 
piecewise-linear model. Whereas the bond extension surpasses the critical extension of unity in the piecewise-linear 
model, and then returns to unity before cracking, the bond extension in the nonlinear model is monotonic, and spends 
a long time near critical extension. A clearer picture of this emerges in Fig. 0t>, where we present the force exerted 
by the bond as a function of time for the nonlinear model. We see that the bond exploits the nonlinearity of the 
force-law to crack in a very slow fashion. 



IV. THE LARGE VELOCITY REGIME 



It is known |5| [l]J that the piecewise-linear solution is inconsistent at large velocity due to the cracking of other 
bonds. This has been extensively studied for the Mode III problem in an infinite square lattice and well as for Mode 
I in a triangular lattice ||. We thus expect that for velocities larger than the critical one at which the inconsistency 
first appears, the solutions of our nonlinear model will diverge significantly from those of the piecewise-linear model. 

To begin, we use our steady-state code to find large velocity solutions of the equations of motion. We present in 
Fig. data for v versus A/Aq at three values of a, together with the analogous results from the piecewise-linear 
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FIG. 3: a) Bond extension along crack surface versus t, for a = 40, 100 and 200 together with the piecewise-linear model 
results, b) Bond force along crack surface versus t, for a = 40, 100 and 200. In both figures, N = 3, rj = 0.25, v — 0.1 and 
dt = .l. 



model. We see that for velocities less than about 0.6c, the curves basically coincide, with the smallest a curve lying 
slightly below the others. However, in the range 0.6c < v < 0.7c, the data exhibit a marked sensitivity to a. The low 
a curve is flat in this region, and as a is increased, monotonicity is destroyed and a local maximum and minimum 
are created. These features are completely absent in the piecewise-linear data, indicating that the additional bond 
breaking, absent by construction in the piecewise-linear solution, is responsible. Clearly, once a local maximum in v 
versus A sets in, as it does for large enough a, the steady-state solutions are unstable for A's past the maximum. We 
address this issue in the next section, where we perform a stability analysis of the steady-state solutions. Nevertheless, 
it is interesting to note that past this band of velocities, the dependence on a is again very mild, but very different 
from the piecewise-linear results. 




FIG. 4: v/c versus A/Ag for N = 8, r\ = 0.1. Data is presented for the cases a— 5, 15, and 100, and well as for the 
piecewise-linear model. 
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V. STABILITY ANALYSIS 

Next, we consider the linear stability of the traveling wave state given by the aforementioned Slepyan ansatz 
U> x y(t) = Uy(t — v / x). Specifically, we assume a solution of the form 

u x , y (t) = u y °\r)+e-^Su y (r) (5) 

where r is the traveling wave coordinate t — x/v and u^ is the previously determined solution. For the stability 
problem, we take Su << 1 and expand to linear order. Note that the last term can be written in the alternate form 
e uvt 5u y (r) with 

Siiy (t) = 5u y {T)e VUJT (6) 

Hence, for stability we must have Re uj < 0, i.e. perturbations must decay in time at a fixed position in the frame 
moving with the crack tip. 

We substitute this assumption into the equation of motion. This leads to the linear problem 

^Su y (r) = £G(f,f) (Su y (T)-6u y ,(T')e^- x '>) 

nn 

+Vk S , S ' (u y °Hr) - t#V)) J: (<Kto - 

where we have defined 

G(g, x ') = /' (uf (r) - uf (r'))+ (<> (r) - uf (r')) ± «> (r) - uf (r')) (8) 

and of course r' = < — x'/f . 

To find the allowed values of w, we proceed as follows. First, we impose the asymptotic boundary conditions that at 
all x outside of our explicit lattice points, Su = 0. We then pick an arbitrary normalization by fixing Su y= ^^(0) = l 
and simultaneously relax the equation of motion at that same point. This procedure, implemented for some guessed 
value of uj, converts our problem into an inhomogeneous, banded system of linear equations for the complex variables 
5u v (t) which can be solved by standard techniques. The missing equation then forms a complex mismatch function 
whose zeroes then determine the allowed values of the eigenvalue uj. 

Two other details should be noted. First, we do not impose any symmetry restrictions on our perturbation with 
respect to reflections about the crack plane. Also, any time we find a root with some specific value of Im uj, an 
equivalent solution exists with 

Im uj — > 2"7r — Im u> 
Su v —>■ Su* (even rows) 

Su v —>■ —Su* (odd rows) (9) 

By even and odd, we mean the parity of the lattice units of distance from the row y = a/3/4. This result follows from 
the fact that the governing equation is real and that the x spacing of nearest neighbor points on the rows separated 
by an even number is ±1, whereas it is ± i for odd rows. A corollary of this is that a root at exactly Im to — tt has 
an eigenfunction which is purely real on even rows and purely imaginary on odd ones. 

One obvious test of our stability analysis arises from the fact that translation invariance guarantees that there is a 
root at to — with eigenfunction 

Actually, our discretiztion for numerical purposes of the continuous variable r breaks this exact invariance and hence 
the mode should be at zero only in the small rfr limit. In the data we present below, the zero mode is found with a 
typical accuracy of 10 -2 , and the eigenfunction to the same accuracy agrees with its expected form. 

In Fig. ||a, we show the basic result that emerges from our analysis, namely that as the driving displacement is 
increased, there is a mode which crosses the stability threshold. In Fig. 5b, we show the associated Imu). At the 
point of instability, Imui is close to but not exactly equal to 7r. Were it equal to 7r, this would be a period-doubling 
instability, as the perturbative displacements at neighboring x values (at the same r) would alternate in sign. More 
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FIG. 5: Eigenvalue for stability problem, with N=8, a — 15 and r\ = .1. The calculation was done with L — 200 and nt=12; 
note that by definition v = l/(ribdt). Re u> is shown in (a), Im ui in (b). 



general values of Im uj signify a Hopf bifurcation which in general has a frequency incommensurate with the original 
frequency associated with moving one lattice spacing. The signature of such a bifurcation should be a slowly oscillating 
crack tip speed, with the oscillation frequency uju p = \Im u) — 7r|. 

In Fig. ^, we show the real part of the eigenfunction for the rows y = ±\/3/4 as a function of the traveling wave 
coordinate r. Note that the perturbation is almost, but not exactly, antisymmetric around the crack plane. The 
eigenvector decays slowly downstream and very rapidly upstream of the crack tip; it is therefore a localized mode 
connected with the tip dynamics diverging from the pure Slepyan form. 

In previous studies of piecewise-linear models (which correspond to the limit a — > oo), it was noted that the 
Slepyan solution breaks down above a critical velocity. There, the criterion for this breakdown was connected to the 
fact that certain bond extensions other than those along the presumed crack line go above the breaking threshold. In 
our current model where the force-law is analytic, there is no such criterion and the steady-state solution found by 
our procedure does not need to be checked for any auxiliary condition. On the other hand, we do observe a linear 
instability. It is therefore of interest to compare these two criteria, i.e. to investigate whether the stability criterion 
in the large a limit goes over smoothly to the inconsistency criterion for the ideally brittle case. To do this, we have 
plotted in Fig. [7] the horizontal bond extension along the crack line for a steady-state solution at a = 50 right at 
the onset of the linear instability. Note that the bond extension is extremely close (roughly of order 1/a) to what 
would be the breaking threshold. Thus, the results of the ideally brittle case as to where nontrivial spatio-temporal 
dynamics of the crack tip sets in are not an artifact of the force-law discontinuity; i.e., smoothing the spring law does 
not alter the basic conclusion. This is good news, as the ideal case has proven rather amenable to analytic techniques 
which enable calculations to be done for much larger N (even for infinite N) than is possible for any direct numerical 
scheme. 

One nagging question posed by our results concerns the fact that the instability is so close to but not exactly at 
the period-doubling point of Im u> = n. We have investigated whether there is any tendency for the mode to get 
pinned at the period-doubling point as parameters are varied, as certain limits are reached etc.. For example, in Fig. 
g we show Im ui at threshold versus a. As far as we have been able to determine, this pinning does not occur and 
the generic bifurcation always leads to an incommensurate oscillating tip state. We will compare this prediction with 
direct numerical simulations of the equations of motion in the next section. 



VI. COMPARISON TO SIMULATION 



In order to test the results of the stability analysis, and also to investigate the nature of the dynamics past the 
instability threshold, we have implemented a direct numerical simulation of the equations of motion. Our system is a 
triangular lattice, of width 2N + 2 whose top and bottom rows are fixed to have displacement ±A. We start with a 
uniformly strained state, except for the first few left hand columns of mass points. The initial displacement of these 
points is taken to vary linearly in x from ±A at x = to the that of the uniformly strained state at x — 10. The 
equations of motion are then stepped forward in time, using a simple Euler scheme. We track the velocity of the crack 
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FIG. 7: Horizontal bond extension for marginally stable solution, with N=4, a = 50 and r\ = .2. The calculation was done 
with L = 200, nb = 12, and dt = .08702. 
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FIG. 8: Variation of Im u) at the marginal stability point with a; data is all for N=4, n — .2. 



by monitoring the bond extensions, and noting when they exceed unity, which we take as our criterion of cracking. 
Note that in this model, cracking is in principle a reversible process and so this criterion is merely a convenient way 
of keeping track of the dynamics and not an intrinsically important threshold. 

As expected, for moderate A the only bonds that "crack" are the diagonal bonds that span the midline. The time 
between cracking events is constant, and, to best compare with our steady-state calculations above, we fix our time 
step dt in each case so that there are 6 time steps between each cracking event; this can then be directly compared to 
our steady-state solutions computed with time bandwidth n& = l/(vdt) =12. We find that the velocities obtained in 
this way reproduce extremely well the velocities calculated by our steady-state code. 

As we increase A above the critical value for instability calculated in the previous section, we indeed find that the 
pattern of steady-state cracking breaks down. Some horizontal bonds on the crack surface begin to crack, and the 
average velocity of the crack falls dramatically, presumably due to the draining of energy by these extra cracking 
events. It is important to note that these broken horizontal bonds are always behind the crack tip, in accord with the 
calculations of the piecewise-linear model. It is not clear whether there is in fact a discontinuous drop in the average 
velocity. The shape of the drop does appear to be insensitive to the dt chosen, so that it in any event does not appear 
to be a numerical artifact of finite dt. 

As a test of our stability calculation, we have computed the Hopf frequency directly from the numerical simulation. 
We fixed A to lie slightly higher than the critical value for the onset of the instability and measured, at the moment 
of the breaking of a northeast- southwest diagonal bond, the extension of the bond. Due to the discreteness of the 
time step, this is somewhat larger than unity. If we had true steady-state propagation, this value would be constant. 
Instead, for this A, we find that the value oscillates about some value, with the magnitude of the oscillation growing in 
time. The growth we associate with the small positive growth rate of the perturbation in the unstable state, and the 
frequency of oscillation with the Hopf frequency. In Fig. ^, we show the data. The period of oscillation is essentially 1, 
(in other words, the time between breakings of this type of bond), which corresponds to a Hopf frequency of essentially 
7r. By numerically fitting a sine wave to the data, we find that the Hopf frequency is approximately 3.272, in very 
good agreement with the stability calculation result shown in Fig. 

As we increase A further, at some point the average velocity starts to increase again. Examining the solution in 
this region in more detail, we show in Fig. 10 the time history of three adjacent points on the crack surface. We see 
that the third trace is identical to the first, and differs significantly from the second. It appears then that the system 
has settled into a new kind of steady-state solution, more complex than the simple Slepyan form. To test this, we 
implemented a new steady-state algorithm, designed to allow for this period-2 type solution. We assumed that in each 
row there are two possible time histories for particles, which alternate. We also did not assume any symmetry across 
the crack surface. Taking care to organize the storage so that the banded structure is maintained, we employed our 
banded Newton's method. Working for the sake of computational convenience at N — 4, the algorithm succeeded in 
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FIG. 9: Bond extension at breaking as a function of position x along the crack surface. Here N = 8, 77 = .1, a=15, A = 2.6215, 
A/A G = 1.2693, di = 0.1053852. The fitted curve is 3.57 ■ 10" 5 sin(3.27166x - 2.85642) exp(0.00368a;). 



obtaining steady-state solutions where the symmetry between the two time histories (as well as up-down symmetry 
across the crack surface) is broken. We traced out the velocity-A relation for this solution. We find that the line 
of solutions is disconnected from the symmetric branch and that in general there are multiple solutions for each A. 
Starting from a solution obtained from using the simulation to obtain an initial guess, we tracked the solution as we 
varied A. The velocity decreased with decreasing A and then turned around and after some wandering headed off to 
large A with the velocity higher than the original branch. Thus it is the slower branch which is physically realized. 
This data is presented in Fig. along with the (time-averaged) velocities measured in simulation. 

It is reasonable to assume that the existence of a period-doubled solution very close to the original symmetric 
branch is associated with the fact that the Hopf bifurcation is nearly period-doubling. To check this notion, we note 
that it was found in the previous section that the Hopf-frequency passed through n for some particular value of a. 
For this value, one expects the period-doubled branch to hit the main branch, corresponding to a higher co-dimension 
bifurcation. To address this question, we have also traced out the bifurcated solution for a = 40, which is very 
slightly below this crossing point. The amazingly baroque solution curve is presented in Fig. |l2|, together with the 
symmetric branch. We see that the bifurcated solution curve turns up to approach the symmetric curve, before 
veering away. Examination of the solutions indicate that the asymmetry is extremely small in this region of close 




FIG. 10: Time trace of displacement of three adjacent points on the upper crack surface. Here N = 8, T) = .1, a=15, A = 2.7, 
A/A G = 1.3073, dt = 0.10858. 
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FIG. 11: Velocity vs. A/Ac for symmetric and period-2 solutions, together with measured (time-averaged) velocity from 
simulation. TV = 4, r\ = 0.1, a — 15. 



approach. Presumably, then, at the crossing value of a the bifurcated solution actually meets the symmetric branch, 
and the bifurcation is perfect at this point. As the meeting is not on the presumably stable branch of the bifurcated 
solution, the bifurcation should prove to be subcritical. A detailed study of the nature of the bifurcation should prove 
mathematically interesting, though not necessarily physically relevant. It should be noted that we have also uncovered 
another, apparently disconnected, branch of period-2 solutions. These other solutions presumably play an important 
role in the exact nature of the bifurcation, though again they do not seem to be important for the physics. 

Increasing A further, the simulation tracks the bifurcated steady-state solution until at some point yet additional 
bonds are broken, and the velocity falls below that of the bifurcated solution. The time-dependence of this state 
becomes more complicated than the simple period-2 structure of the bifurcated steady-state. It is important to point 
out that the additional bonds are also always located behind the crack tip, which remains on the midline. Thus the 
additional bond breaking is reminiscent of side- branching in dendritic growth [fjjsf ; i.e., a phenomenon induced by 
the tip and which is left behind by the growing crack. As opposed to dendritic growth, in Mode III fracture the 
sidebranching is generated intrinsically by the tip, and is not noise-induced |rij] . We see, however, no tendency for 
the tip itself to leave the midline. This can be seen in Fig. [l^, where the broken bonds are portrayed for two different 
values of A. Also, the picture shows that as the driving is increased the size of the sidebranches increases, but they 




FIG. 12: Velocity vs. A/Ag for symmetric and period-2 solutions. TV = 4, r\ = 0.1, a = 40. 
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FIG. 13: Broken bonds for N = 20, 77 = .1, a = 15. (Top) A=5.5, A/A G =1.715; (Bottom) A=6.5, A/A G =2.207 



are always microscopic, growing to a length of about 8 at the larger A. In this latter picture, the crack is moving 
at an average velocity of v/c = .789, which is quite fast. Furthermore, the sidebranches are not only short, but the 
sidebranching period is also on the lattice scale. Thus, the claim of |ll[ that these sidebranches are related to the the 
experimentally seen micro-branching appears to us doubtful, especially in light of our preliminary work on Mode I 
cracking fisj| , where the tip dynamics is very different. This issue is worthy of further exploration. 
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VII. APPENDIX - SOLVING THE ALMOST BANDED PROBLEM 



In this appendix, we outline how to efficiently solve the linear problem associated with the Newton's method 
treatment of the nonlinear steady-state problem. It is important to set up the problem correctly to achieve an almost- 
banded structure. We order the variables with the y index running more quickly, since we want N -C 2L + 1, the 
length of the system in the x direction. We order the equations similarly, replacing the equation of motion for the 
site x = 0, y = Vo/4 (where x runs from — L to L) by the constraint equation u = uO, where the constant uO can be 
chosen arbitrarily. The equation of motion for the site (0, a/3/4) is taken to be the last equation. We also note that 
equations of motion have to be taken for the columns — L — 1 < x < L — 1, since the number of linear modes which 
diverge as x — > —00 is one greater than the number which diverge for large positive x (see || [7J) due to the presence 
of the third derivative operator. 



The problem then has the structure 



M b V\ ( SU\ _ ( A 

w T J Ua J - I b 



(11) 



Here, Mb is a banded square matrix of size N(2L+1) with lower bandwidth N(nb + 2) and upper bandwidth N(nb + 1) 
(with rib an integer such that the velocity is v = l/(ribdt)). The replacement of the equation for the distinguished 
site (0, %/3/4) is necessary to ensure that M is not singular. The length N(2L + 1) vector V contains the derivatives 
of the equations of motion with respect to A and the vector W contains the equation of motion for the distinguished 
site (which does not depend on A). The shift in the displacements SU, and the souurce A are also vectors of length 
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N(2L + 1), and 5A and b are numbers. We can explicitly solve this system in terms of M. h 1 , which is easy to compute 
due to its banded structure. We find 

W T MZ 1 A-b 

<5A = =-2 ; 

W T M^V 

SU = M^(A - (SA)V) (12) 
Hence, our entire system can be solved with no more effort than would be required for a fully banded problem. 
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